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Abstract 

We numerically solve the propagation of a shock wave in a chain of elastic 
beads with no restoring forces under traction (no-tension elasticity). We find 
a sequence of peaks of decreasing amplitude and velocity. Analyzing the main 
peak at different times we confirm a recently proposed scaling law for its decay. 
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1 Introduction 



The study of granular matter is of interest for fundamental physics as well as for 
technological applications. The intensive studies of liquids and solids gave physicists 
powerful tools for investigating these states of matter. These tools are however very 
difficult to apply to granular matter jj| . 

Few works have been done on sound propagation in real granular media. Exper- 
iments have been performed by Liu and Nagel |^|, ||, ^] who investigated the prop- 
agation of low amplitude vibrations in a box of spherical beads. It was concluded 
that non-linearity and disorder make wave propagation behave in very unexpected 
ways: the fragile structure of contacts between grains make them sensitive to re- 
arrangements changing dramatically the amplitude response of the receptor to the 
source. Besides, disorder gives rise to important interference effects that can lead 
to localization || . Numerical simulations have been made by Melin Q who studied 
the depth dependence of sound speed in granular media. The results obtained were 
different from the power law behavior between sound speed and pressure predicted 
by Goddard § in an effective medium calculation. 

Three dimensional models seem to be very difficult to investigate directly and 
we need to begin with simplified models in order to isolate specific features of real 
granular media. In this paper we try to understand the effect of non-linearity on 
wave propagation in one dimension. A previous study was done by Nesterenko Q 
on a chain of spherical beads obeying the elastic Hertz law of contact. It was shown 
both analytically and numerically that the chain submitted to compression pulses 
involves solitary wave propagation. This has been confirmed by the experiments of 
Lazaridi and Nesterenko |9| and of Coste, Falcon and Fauve |l(| . 

Here, we want to study the case of a chain of beads initially in contact and 
submitted to a shock perturbation. The beads interact via an elastic contact law 
only when they are compressed. The problem is therefore highly non-linear because 
as soon as there are broken contacts the chain behaves as an ensemble of independent 
elastic systems. We can look at this problem even in terms of a linear chain of 
points connected by springs with completely asymmetric elastic constants: k in 
compression and zero in extension. Here we must precise that the Hertz law which 
has been studied until now correspond to perfect spherical beads. In a real granular 
medium the shape of the beads at contacts can be far from spherical. In fact beads 
interactions can be modelized with a force law exponent varying from one to two. 
Our linear compression force model belongs to this range. This one-dimensional 
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system is clearly far from being realistic but it exhibits some features which can help 
understanding of the more general problem. In particular we try to elucidate how 
the front wave propagates and how its lost energy can contribute to the formation 
of several other secondary waves. The outline of the paper is as follows: in section 
2 we define our model and we present the setup for the numerical simulations. 
Section 3 is devoted to the description of the results concerning the phenomenology 
of the system. In section 4 we discuss in detail the mechanism of front formation. 
Discussions and conclusion will be given in section 5. 

2 Model definition 

Our model is composed of N spherical beads of radius R and mass m. We define 
the force between two neighboring beads as varying linearly under compression and 
being equal to zero when the beads are not in contact. It can be written as: 

f = kse(d) (i) 

where k is the spring constant, 2R — S the distance between the centers of the two 
neighboring beads and 9 the Heaviside function. Mechanically such systems are 
known as "no-tension elastic" fj^ ]. Initially all the beads are just touching each 
other without exerting any forces on each other except for the first bead which 
penetrates the second one by So . We then let the system evolve with the left 
end fixed and the right end free. An experimental realization of the model can 
easily be done as it can be seen from the Fig. ([!]). This one dimensional array of 
non-connected linear springs is the simplest model for no-tension elasticity. 
The equations of motion are the following: 

mui = k5i-i,i9(6i-i,i) - k8 i>i+1 9{8 i4+1 ) 0<i<N-l (2) 

where Ui is the displacement of bead number i. One can write: 

<5i_ M = 2R - (xi - x^i) (3) 
= 2R - (ui + x i>0 - + ar»-i,o)) ( 4 ) 
= m-i - Ui (5) 

where Xi is the position of bead number i and Xi t o — (2i + 1)R. 
Thus, 

mil, = k(v,i-i - Ui)6{Ui-\ - m) - k(ui - u i+1 )9(ui - u l+1 ) < i < N - 1 (6) 
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with the boundary conditions given by 

Un = Sn for t > 

(7) 

Ui = for i > and t = 

These boundary conditions correspond to the study of the propagation of a shock 
in the chain. From now onwards we will denote with forward and backward di- 
rection the direction of increasing and decreasing i, respectively. For the nu- 
merical implementation of the analysis we have chosen N = 1500, R = bmm, 
p = 1.9 x 10 3 fcg.m -3 , k = 10 7 N.m^ 1 and 5q — 0.5mm. Several algorithms have 
been used to solve the system of equations (6). First, we implemented the Verlet 
scheme which was not able to give a good enough precision. Much better accuracy 
was obtained with a 5th order Gear predictor-corrector method being even more 
precise than the 4th order Runge-Kutta scheme. The test of the accuracy was 
obtained by monitoring the energy conservation of the system. During the compu- 
tation, energy conservation was satisfied with a relative error lower than 10 -6 with 
a time step At = 10 _8 s. It is worth to stress that during the total time of evolution 
the perturbation never reaches the end of the chain. 



3 Front phenomenology and comparison with the 
harmonic case 

In order to follow the evolution of the perturbation we have monitored the evolution 
of the displacements Ui and of the velocities Vi of the beads versus the bead number 
i as a function of time. A snapshot of the solutions at a given time is shown in 
Fig.(|2|) and Fig.(|J). The bead displacements are characterized by a front wave 
followed by inclined plateaus. The plateaus consist of an almost smooth curve 
with jumps. It is worth to note how, in the plot of the displacements, intervals 
with positive slopes typically correspond to regions of detached particles while the 
negative slopes correspond to particles in contact with each other. The snapshot for 
the velocities exhibits a structure of peaks with decreasing amplitude. These peaks 
correspond to waves propagating with monotonically decreasing velocities. Each 
wave is, at its turn, composed by several particles which are all in contact with each 
other and moving in the same forward direction. In between the peaks particles are 
not in contact and propagate in the backward direction with decreasing velocities. 
In Fig. (||) we also see the appearance of noise. This is a real numerical noise which 
is reduced but never disappears when we decrease the time step resolution. In fact, 
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decreasing of one order of magnitude the resolution allows us to observe another 
peak and so on and so forth. In the ideal case of infinite resolution one would be 
able to observe an infinite series of peaks behind the front. As we shall see later this 
fact has no consequences on our conclusions since we will discuss the region where 
no noise occurs and results can be extrapolated to the noisy region. Coming back 
to the analysis of Fig.(^) we can say that periodically the front looses energy when 
its last particle is detached, i.e looses contact. This is the only way in which the 
energy stored in the front wave can decrease. In this way behind the front one has 
a chain of particles moving backward with decreasing velocity. These particles can 
eventually contribute to the formation of peaks whose nature is different from that 
of the front. We will come back to this point later when we will discuss the actual 
mechanism for the front formation. 

It can be easily shown that the front wave propagates nearly at the sound veloc- 
ity. In order to do that we calculate analytically the speed of sound in the medium 
of the same density of our beads and we compare the results with the speed observed 
numerically for several values of the elastic constant k: 10 8 N.m , 10 7 N.m^ 1 and 

io 6 N. m - 1 . 

The calculation of the speed of sound is made as follows: if we consider the 
beads as springs of spring constant k and length 2R we get a dispersion relation 
that is: 



where q is the wave vector. The speed of sound being defined by c s — lim g _,o ^, we 
get 



Replacing all the parameters with their numerical values we have c s = 3170. 47m. s 
when k = 10 s N.m' 1 , c s = 1002. 59m. s _1 when k = l^N.m^ 1 and c s = 317.05m.s~ 1 
when k = lQ G N.m~ l . These speeds are systematically slightly greater than the ones 
computed numerically which gave us 3160m. s -1 , 1000?ti.s _1 and 316m. s -1 respec- 
tively. These numerical values have been obtained by counting the number of beads 
separating the positions of the maxima of velocity in the front at two different times. 
We cannot thus pretend an exact agreement with theoretical values. 

In order to have a better understanding of our results we compare our system 
with the harmonic chain consisting in a series of springs with spring constant k. 
We integrated the harmonic chain keeping the same initial conditions. After some 
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algebra we get the chain eigenfrequencies: 

^ 2 = 4- cos 2 — ^— 1 < Z < JV-1 (10) 
m 2JS — 1 

We then obtain for the displacements: 

i u n (t)= S Q + (-l) n J2f l =1 1 C l sin cos u t t 1 < n < N - 1 
[ u (i) = 5o 

where the C; are constants depending on the initial condition. 

This relation is shown in Fig. (|J) where we recognize the well-known oscillations 
after the front due to the discretisation as already discussed by Gibbs and being now 
what we call the "Gibbs phenomenon". We also plotted in Fig.([|) what happens 
when one tries to go continuously from the harmonic case to the no-tension case 
by decreasing the spring constant of the springs when they are under traction and 
keeping them equal to k under compression. Each curve corresponds to a different 
value of the ratio of the spring constant under traction over k. We clearly observe a 
continuous transition from the harmonic regime to the no-tension one as we decrease 
the ratio of the spring constants. 



4 Wave formation 

In the previous section we have made some observations about the displacement 
and velocity profiles for a fixed time of propagation. Let us now consider the time 
evolution in order to understand how the peaks are created and how they propagate. 

4.0.1 Scaling of the front shape 

We define the amplitude of the front wave as the maximum velocity of the beads 
belonging to it. The curve representing this amplitude as a function of time is 
shown on Fig.(^). One can see that the amplitude decreases and oscillates. These 
oscillations are due to the discretisation, i.e to the fact that the bead having the 
maximum velocity keeps it during a finite time corresponding in Fig.(^|) to the width 
of the oscillation. Initially the fastest bead has a low velocity that increases towards 
a maximum with time. Then this value decreases to a value lower than the initial 
value. This behavior continues with the right neighbor of the previous bead which 
then has the new maximum speed and keeps it during a finite time with the same 
evolution as before. On a log-log scale we find that the envelope of the maxima is 
well fitted by a power law A(t) oc t~ a with a ~ 0.17, consistent with a t~i behavior. 
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This value of the exponent seems to be universal since it is robust when changing 
the value of the elasticity constant k in the range 10 6 — 10 8 and Sq in the range 
0.3 — 0.7mm. 

We also define the width of the front wave by counting the number of beads 
with velocities greater than zero belonging to the front. For this quantity we find 
a scaling law L(t) cx t@ with ft ~ 0.32 (see Fig.(^)). In this case the behavior is 
consistent with ti and is very robust with the same universal character as for the 
scaling of A(t). The exponents a and ft agree with the calculation made by J. Hinch 
fPUl using the conservation of kinetic energy of the front wave (a' = ^,ft' = |). 

The use of the two previous power laws enables us to find a scaling relation 
for the front. By plotting ^r=Q versus x ' t/ ^ :ft where vj is an adjustable parameter 
with the dimension of a speed, we found that the curves corresponding to different 
values of t collapse on a single one as shown on Fig.(0). The collapse happens for 
a value of Vf = 1001. 5m. s -1 which is very close to the sound speed 1002. 5m. s -1 
(k = 10 7 N.m' 1 ). It is worth to note that v* is the speed measured at the center 
of mass of the front and not at its maximum. The best collapse is obtained for the 
values of a = 0.17 and ft' = |. The exponent ft changed to ft' in agreement with 
the theoretical results whereas a is not exactly equal to a' since its calculation is 
based on kinetic energy conservation. In the present case energy is dissipated from 
the front wave by the effect of detachment of particles which induce a correction to 
a' . We can therefore write the following scaling- law for the elastic front: 



where / is a scaling function that has also been calculated |ll|] and Vf = 1001. 5m. s . 

This scaling relation allows us to describe the self-similarity of the front wave 
shape. Another interesting feature coming from the scaling relation is that the 
velocity of the front is defined at the center of mass of the beads belonging to it. 
The other points on the front shape move with different velocities. Each point 
behind the center of mass moves with a velocity smaller than the sound speed 
whereas the points on the part preceding the center move with a larger velocity. 
Thus the front wave is divided into a subsonic and a supersonic part. This kind of 
wave is very different from the solitonic ones discovered by Nesterenko Q where it 
was found that the wave travels with a constant shape and a constant speed for all 
its points. 




(12) 
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4.0.2 Emergence of the secondary peaks 

Let us now focus our attention on the secondary peaks propagating behind the elas- 
tic front. The particles detached from the front move backwards and their energy 
does the same (the system is globally conservative). This energy propagating back- 
wards is then transferred from particle to particle until it reaches the left boundary 
where it is reflected and starts to propagate forward. This is the beginning of the 
propagation of a peak. The process then continues by means of detachments of 
particles from this peak contributing in this way to the formation of new peaks. 

In order to analyze in detail the evolution of the peaks we looked at the spatio- 
temporal structure of the contacts of the first one hundred beads of the chain. This 
can be seen in Fig. (||) : a grey dot means that there is a contact between two beads 
whereas a white dot means that they are not in contact. The computation of the 
contact distribution of the beads is done every 10 s up to a propagation time 
of 2.5 x 10 _3 s. At lCP 5 s almost all the beads are under compression. Then this 
compression chain moves to the right letting empty spaces between the first beads of 
the chain. This corresponds to the formation of the front wave. As the compression 
chain travels (the large grey triangle on the right lower corner of the plot) one can 
observe grey dots appearing at the beginning of the chain. This means that a new 
compression wave is being created. This wave is much thinner than the front wave 
and is represented on the plot by an almost straight curve with a slope slightly 
above the slope of the front showing that this first peak travels slower than the 
front. It can be seen that the other peaks are also formed at the begining of the 
chain all making an almost straight grey line with a slowly increasing slope which 
means a decreasing speed. It is easy to understand that the existence of these peaks 
is due to reflections on the first bead which is maintained fixed playing the role of a 
reflecting wall. Hence without noise there should be more than five peaks in fig.(|J): 
all the peaks behind these five ones have been destroyed by a numerical noise and 
one should see a succession of peaks with decreasing amplitude until the beginning 
of the chain. One more important point emerging from the spatio-temporal pattern 
is that beads which do not belong to a peak are not under compression and are thus 
in a ballistic regime as one can see from the white spaces separating the peaks in 
Fig.®. 
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5 Conclusion 



We have studied the response of a chain of beads to a small displacement when 
no static force is applied to its ends. We observe the propagation of waves very 
different from the solitons which can be seen in a hertzian chain. 

The non-linearity of the problem lies in the contact law used: a step function 
which indicates that there is no force when two beads are not in contact. 

By keeping the left end of the chain fixed we observed an interesting bounc- 
ing phenomenon. In addition to the elastic front, secondary peaks of decreasing 
amplitude are generated at the left end of the chain. These peaks correspond to 
several beads in compression traveling in the forward direction whereas the beads 
in between the peaks travel in the backward direction without touching each other. 
We have found an interesting scaling law for the velocity of the beads belonging to 
the elastic front. This law which has also been derived analytically by J.Hinch |llj 
seems to be universal since it does not depend on the parameters of the problem. 

Starting form these results it would be interesting to look at a chain submitted to 
other types of perturbations. Besides, it is important to look at different networks 
(two or three dimensional) in order to understand the effect of the network structure 
on wave propagation. 
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Figure 1: Simple experimental realization of the model. 
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Figure 2: Displacements of the beads at propagation times equal to 6 x 10 _3 s, 10~ 2 .s 
and 1.4 x 10~ 2 s. 
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Figure 3: Velocities of the beads after 1.4 x 10 2 s. 
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Figure 4: Curves of displacements for different values for the ratio of the elasticity 
constants. The ratio goes from one (harmonic chain) to zero (nonlinear chain). The 
propagation time is 10~ 2 s. 
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Figure 5: Double-logarithmic plot of the amplitude of the front wave versus time. 
After a fast increase it decreases and oscillates. The dashed line corresponds to the 
function f(x) oc x~ a which fits the envelope of the maxima of the amplitude curve. 
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Figure 6: Double- logarithmic plot of the width of the front wave versus time. The 
dashed line corresponds to the function f(x) oc x 13 showing an increase of the width 
with a power law. 



14 




Figure 7: Collapse of the front wave shape to a single curve verifying the scaling 
law for the bead velocities in the front wave. We show four curves corresponding 
to t = 0.3 x l(T 2 s,0.6 x lCT 2 s,0.9 x 10" 2 s,1.2 x lCT 2 s. a = 0.17, P' = § and 
v f = 1001. 5m.fi- 1 . 
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Figure 8: Spatio-temporal evolution of beads contacts. The first one hundred con- 
tact points of the chain are on the horizontal axis and time is on the vertical axis. 
The distribution of contacts is computed every 10~ 5 s. Here we started at 10~ 5 .s 
till 2.5 x 10~ 3 s. A white square means that there is no contact between two beads 
whereas a grey one means that the beads are in contact. 
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